% Specify x range and number of points
x0 = -3;
x1 =  3;
Nx = 301;

% Specify y range and number of points
y0 = -3;
y1 =  3;
Ny = 301;

% Construct mesh
xv    = linspace(x0,x1,Nx);
yv    = linspace(y0,y1,Ny);
[x,y] = meshgrid(xv,yv);

hold on;

% Calculate z
z = x + 1i*y;
% 1nd order Runge-Kutta growth factor

% Calculate magnitude of g
g1 = 1+z;
gmag1 = abs(g1);
% Plot contours of gmag1
contour(x,y,gmag1,[1 1],'k-','linewidth',2);

% 2nd order Runge-Kutta growth factor
g2 = 1 + z + 0.5*z.^2;
gmag2 = abs(g2);
contour(x,y,gmag2,[1 1],'r-','linewidth',2);

% 3nd order Runge-Kutta growth factor
g3 = 1+z+0.5*z.^2+(1/6)*z.^3;
gmag3 = abs(g3);
contour(x,y,gmag3,[1 1],'b-','linewidth',2);

% 4nd order Runge-Kutta growth factor
g4 = 1+z+(1/2)*z.^2+(1/6)*z.^3+(1/24)*z.^4;
gmag4 = abs(g4);
contour(x,y,gmag4,[1 1],'g-','linewidth',2);

legend({'RK1','RK2','RK3','RK4'},Location="best");
legend off;
axis([x0,x1,y0,y1]);
axis('square');
xlabel('Re(\lambda\Delta t)');
ylabel('Im(\lambda\Delta t)');
title('Stability regions of Runge-Kutta');
grid on;grid minor;box on;